Method for fast Fourier transformation for nuclear magnetic resonance signals arranged in a matrix having an arbitrary matrix size

ABSTRACT

A Fourier transformation in a measured data matrix (MD) with the discrete measured data (x j ,k) is implemented for all k=0 . . . N-1, respectively in the following steps, whereby it is defined: w=e -2i π/N 
     a) The measured signals (x j ,k) are acquired with an rf-phase modulation of the nuclear magnetic resonance signals according to a chirp function w j .spsp.2 /2 , 
     b) The arising data are convoluted with a second chirp (w - (j-k).spsp.2 /2 ) and entered into an auxiliary data matrix (HD), 
     c) The image data matrix is acquired from the auxiliary data matrix (HD) by Fourier transformation in at least one further direction.

FIELD OF THE INVENTION

The invention is directed to a method for image acquisition frommeasured data of a nuclear magnetic resonance experiment that areentered into a measured data matrix with discrete k-space points,whereby a image matrix with image data is acquired from the measureddata matrix by Fourier transformation, and whereby an image is displayedon the basis of the image data matrix.

DESCRIPTION OF THE PRIOR ART

A method of the above general type is disclosed, for example, in U.S.Pat. No. 5,168,226.

In traditional methods, the Fourier transformation required for imageproduction is implemented according to what is referred to as an FFT(Fast Fourier Transform) algorithm, as described, for example, by Alfredv. Aho et al. in "The Design and Analysis of Computer Algorithms",Addison-Wesley Publishing Company, pages 257 through 264. This algorithmis very efficient since it requires only a fraction of the calculatingoperations compared to discrete Fourier transformation, however, it canonly be applied when the length of the dataset to which the algorithm isto be applied represents a power of a whole number. The FFT method isusually applied to datasets having a length of 2^(n).

The limitation to specific dataset lengths, however, is undesirable forcertain applications. For three-dimensional datasets, for example, theresolution within a slice should thus be capable of being arbitrarilyset. This can in fact be achieved without further difficulty with adiscrete Fourier transformation, however, the number of requiredcalculating operations and thus the required calculating time, increasesdrastically. For example, given a dataset length of 512 more than onehundred times as many calculating operations are required for thediscrete Fourier transformation than given the FFT algorithm.

Chirp-z transformation as disclosed in European Application 0 453 102 inits application to MR imaging is a method for image calculation that isindependent of the dataset length. This method is noticeably faster thandiscrete Fourier transformation.

Utilizing a "Chirp-Z transformation" for the realization of a discreteFourier transformation is known for certain applications such as, forexample, "charge coupled devices" from the reference "Henry Nussbaumer,Fast Fourier Transform and Convolution Algorithms", Springer-Verlag,1982, pages 112-115.

The technique of "Chirp-z Transformation" is explained in general fordigital signal processing in the book, "Digitale Signalverarbeitung",Vol. I, H. W. Schussler, Springer Verlag 1988, pp. 70-72.

The reference, "IEEE Transactions on Medical Imaging", Vol. 9, No. 2,Jun. 2, 1990, pp. 190-201, describes the application of a chirp-ztransformation for chemical shift imaging with magnetic resonanceapparatus. This transformation is thereby utilized as substitute for theFourier transformation for acquiring the spectrum, and an improved,spectral resolution is thus achieved. The Fourier transformation isutilized for the topical resolution, as usual.

The reference "Electrical Design News", Vol. 34, 1989, pages 161-170describes the application of chirp-z transformation to the calculationof frequency spectra.

SUMMARY OF THE INVENTION

An object of the present invention is to implement a method for imagereconstruction in a nuclear magnetic resonance tomography apparatus suchthat measured data matrices of an arbitrary size can be processed givenlow calculating time.

The above object is achieved in accordance with the principles of thepresent invention in a method for image acquisition from nuclearmagnetic resonance measured data (x_(j),k) which are entered into ameasured data matrix with discrete k-space points (j,k), with a Fouriertransformation being implemented in a first direction for all k=0 . . .N-1 (wherein N is the number of rows in the matrix) and wherein themeasured signals (x_(j),k) are acquired with an RF-phase modulation ofthe nuclear magnetic resonance signals according to a chirp functionw^(k).spsp.2^(/2), wherein the resulting data are convoluted with asecond chirp function w⁻(j-k).spsp.2^(/2), and are entered into anauxiliary data matrix, wherein an image data matrix containing the imagedata (x_(j),k) is acquired from the auxiliary data matrix by Fouriertransformation in at least one further direction, wherein an image isthen displayed on the basis of the image data matrix, and whereinw=e⁻²π/N.

DESCRIPTION OF THE DRAWINGS

FIGS. 1-5 illustrate a conventional spin-warp pulse sequence forexplaining the problem addressed by the inventive method.

FIG. 6 schematically illustrates the arrangement of a measured datamatrix generated according to the conventional pulse sequence shown inFIGS. 1-5, and used in the inventive method.

FIG. 7 is a block diagram illustrating the basic steps of the inventivemethod.

FIG. 8 is a chart illustrating the convolution steps in accordance withthe principles of the present invention.

FIG. 9 is a block diagram showing radio-frequency modulation inaccordance with the principles of the present invention.

FIG. 10 is a table showing the dependency of the required calculatingsteps on the Fourier transformation method which is used.

DESCRIPTION OF THE PREFERRED EMBODIMENTS

The convention spin-warp sequence of FIGS. 1 through 5 is merelyintended to serve for explaining the problem addressed by the inventionmethod. It must be emphasized that the inventive method can be appliedfor all pulse sequences wherein the information about the spatial originof the signal contributions required for image generation are coded inphase factors of the nuclear magnetic resonance signal. Given theexample of a pulse sequence shown in FIGS. 1 through 5, afrequency-selective radio-frequency pulse RF is first emitted in underthe influence of a slice selection gradient G_(S). Nuclear spins areexcited in only one slice of the examination subject under the influenceof the slice selection gradient G_(S). Subsequently, the dephasingcaused by the positive sub-pulse of the slice selection gradient G_(S)is in turn cancelled by a negative sub-pulse G_(S) ⁻. Further, aphase-coding gradient G_(P) is emitted. Finally, a negative readoutgradient G_(R) ⁻ is also activated.

Only a positive readout gradient G_(R) ⁺ is activated during thefollowing readout phase. The arising echo signal is sampled M times andthe M measured values acquired in this way are entered into a row of ameasured data matrix MD according to FIG. 6.

The illustrated pulse sequence is repeated N times with different valuesof the phase-coding gradient G_(P), so that a measured matrix with Nrows is obtained overall. Usually, the phase-coding gradient is therebyswitched in equal steps from the highest positive to the highestnegative value or vice versa from pulse sequence to pulse sequence. Themeasured data matrix MD can be considered as a measured space, as ameasured data plane in the two-dimensional case of the exemplaryembodiment. This measured data space is referred to as "k-space" innuclear magnetic resonance tomography.

The information about the spatial origin of the signal contributionsrequired for the image generation is coded in the phase factors, wherebythe relationship between the locus space and the k-space mathematicallyexists via a Fourier transformation. Valid is:

    S(k.sub.x,k.sub.y)=∫∫ρ(x,y)e.sup.i(k.sbsp.x.sup.x+k.sbsp.y.sup.y) dxdy

The following definitions thereby apply: ##EQU1##

Given the measured data matrix MD shown in FIG. 6, each row correspondsto a single nuclear magnetic resonance signal S. The sampling in thek-space ensues in successive rows given step-by-step switching of thephase-coding gradient G_(PC). At the beginning of each sub-sequence, aphase-coding gradient G_(P) whose gradient amplitude rises continuouslyin steps from sub-sequence to sub-sequence is respectively activatedbefore the first nuclear magnetic resonance signal S₁. When, forexample, each nuclear magnetic resonance signal is sampled with 256measured points and 256 phase-coding steps are implemented, then a rawdata matrix with 256 rows and 256 columns is obtained, i.e. 256×256measured values in the k-space. The analog measured values obtainedgiven the pulse sequence of FIGS. 1 through 5 are thus digitized onto agrid in the k-space. In the illustrated two-dimensional case, a rownumber j and a column number k is allocated to each digitized measuredsignal, referenced x_(j),k below. In order to obtain a image data matrixBD from the measured data matrix, the measured data matrix isFourier-transformed both in row direction as well as in columndirection. If a discrete Fourier transformation were thereby applied,then the following complex sums would have to be calculated for all k=0. . . N-1 for, for example, the Fourier transformation in columndirection. ##EQU2## whereby X is an image data point and w=e^(-2i)π/N.

The same is true of the Fourier transformation in row direction.However, a considerable calculating outlay is thus involved.

The Fourier transformation can be implemented significantly moreelegantly with the FFT (Fast Fourier Transform) algorithm. Given amatrix size of 512×512 measured points, the number of requiredcalculating operations sinks to one one-hundredth. As already initiallyexplained, the FFT algorithm, however, only functions for arrays whoselength is a power of a number, whereby powers of 2 are normallyemployed. This limitation of the matrix size to specific numericalvalues, however, is disturbing in some instances.

In order to avoid this disadvantage, the inventive method makes use ofthe principle of discrete Fourier transformation, but achieves asolution wherein significantly fewer calculating steps are required thanin the above-recited calculation of the complex sums. Namely, what isreferred as a chirp-Z transformation algorithm (CZT) is applied. Thefollowing algebraic identity is utilized in this transformationalgorithm: ##EQU3##

As shown in the block circuit diagram of FIG. 7, the CZT algorithm isrealized in the following steps:

Each digitized measured signal x is multilied by a chirpW^(j).spsp.2^(/2) (as shall be presented later, this effect is achievedin that the measured signals are acquired with a corresponding phasemodulation).

A convolution ensues with a second chirp w⁻(j-k).spsp.2^(/2).

A multiplication ensues with a further chirp w^(k).spsp.2^(/2).

The most complicated operation is the convolution, whereby an array of Nmeasured values X_(j),k must be convoluted with an array of N datapoints of the chirp w⁻(j-k).spsp.2^(/2) for each column k of themeasured data matrix MD in the illustrated example. The steps therebyimplemented for a column of the measured data matrix are shown in FIG.8. First the two arrays for x_(j),k and w⁻(j-k ).spsp.2^(/2) are broughtto at least double the length by filling with zeroes. For clarity, theillustrated example proceeds on the basis of only four data points thateach are brought to an array length of eight data points by filling withfour zeroes. An FFT algorithm is applied to both arrays. The twoFourier-transformed arrays are multiplied point-by-point. An array isthus obtained to which an inverse FFT algorithm is applied. As a result,one finally again obtains an array with eight data points as result ofthe convolution. Only four of these data points are further-processed;the rest are discarded.

The filling with zeroes in the initial arrays has the purpose ofavoiding over-convolutions in the application of the FFT algorithm. Whenthe lowest value in the illustrated example is employed as zero pointfor the FFT algorithm, then the four lowest values in the initial arrayare further-processed.

Data that correspond to the result of a discrete Fourier transformationare obtained as result. It is thereby important that no demands are madeof the length of the measured data sets in any of the illustrated steps,so that arbitrary matrix sizes can be processed, as in classic, discreteFourier transformation.

With a given matrix length, the multiplication coefficientsw^(j).spsp.2^(/2) and w^(k).spsp.2^(/2) as well as the Fourier transformof the chirp w⁻(j-k).spsp.2^(/2) represent fixed quantities that,expediently, are calculated once before the chirp transformation andstored in tables for later use. Given the illustrated chirp-Ztransformation, one manages with significantly fewer calculatingoperations, measured in floating decimal point additions andmultiplications, than in classic, discrete Fourier transformation. Thecorresponding values are shown compared in the table of FIG. 10 fordifferent matrix lengths for the fast Fourier transformation (FFT), thechirp-Z transformation (CZT) and the discrete Fourier transformation.One can see that, given larger array lengths, the CZT operation isslower by approximately the factor 5 than the FFT operation involvingthe limitation with respect to array length. However, CZT is faster bythe factor 20 compared to the DFT operation that does not involve theselimitations.

Since FFT continues to be by far the fastest method, this will still beemployed in those directions where it remains expedient on the basis ofthe array length and the use of CZT will be limited to those cases thatcannot be processed with FFT. In three-dimensional imaging, for example,it could be advantageous to apply CZT only in the orthogonal directionto a slice where arbitrary thicknesses or resolutions of the slice aredesired. A conventional matrix size of, for example, 512×512 points, bycontrast, can be used in the other directions and, thus, the faster FFTtransformation can continue to be used.

There are various possibilities of shortening the illustrated QZT methodeven further. For example, only a phase rotation is implemented with themultiplication by the chirp w^(k).spsp.2^(/2) One can thus forego thismultiplication step when the phase of the complex signals obtained isnot evaluated but only the amount thereof is calculated.

According to the invention, the first multiplication step, as shown inFIG. 9, is achieved by radio-frequency phase modulation of the nuclearspin signal with the function w^(j).spsp.2^(/2). This, namely,corresponds to a multiplication with the factor w^(j).spsp.2^(/2).

Although modifications and changes may be suggested by those of ordinaryskill in the art, it is the intention of the inventor to embody withinthe patent warranted hereon all changes and modifications as reasonablyand properly come within the scope of his contribution of the art.

What is claimed is:
 1. In a method for image acquisition from nuclearmagnetic resonance signals entered into a measured data matrix having Nrows with discrete k-space points (j,k), with a Fourier transformationbeing implemented in a first direction for all k=0 . . . N-1, theimprovement comprising the steps of:(a) receiving nuclear magneticresonance signals and acquiring a data set of measured signals (x_(j),k)by RF-phase modulating said received nuclear magnetic resonance signalsin a modulator according to a first chirp w^(j).spsp.2^(/2), whereinw=e^(-2i)π/N ; (b) convoluting said data set with a second chirpw⁻(j-k).spsp.2^(/2) to obtain a convoluted data set and entering saidconvoluted data set into an auxiliary data matrix; (c) generating animage data matrix from the auxiliary data matrix by Fouriertransformation thereof in at least one further direction; and (d)displaying an image obtained from the image data matrix.
 2. A method asclaimed in claim 1 comprising the additional step of multiplying theconvoluted data in the auxiliary data matrix for all k=0. . . N-1 with athird chirp w^(k).spsp.2^(/2).
 3. A method as claimed in claim 2comprising the additional step of storing values for the third chirp ina table.
 4. A method as claimed in claim 1 wherein step (c) comprisesgenerating the image data matrix by magnitude formation from a complex,completely Fourier-transformed data set.
 5. A method as claimed in claim1 wherein the step of convoluting with the second chirpcomprises:filling the data set acquired in step (a) with zeroes toobtain a first data array of at least twice a length of said data setand filling said second chirp with zeroes to obtain a second data arrayof at least twice a length of said second chirp; fast Fouriertransforming each of said first and second data arrays to obtain firstand second fast Fourier-transformed arrays; multiplying said first andsecond fast Fourier transformed arrays point-by-point to obtain aproduct and entering the product into a product array; and acquiringsaid convoluted data from the product array by inverse fast Fouriertransformation of said product array.
 6. A method as claimed in claim 1comprising the additional step of storing values for the first chirp ina table before executing steps (a) and (b).
 7. A method as claimed inclaim 1 comprising the additional step of storing a fastFourier-transformed second chirp in a table before executing steps (a)and (b).
 8. A method as claimed in claim 1 wherein a Fouriertransformation in steps (a) and (b) is implemented only for directionsof said measured data matrix that do not comprise a data set lengthsuitable for a fast Fourier transformation, and implementing fastFourier transformations for other directions.